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We present a suite of Mathematica-based computer-algebra packages, termed "Kranc", which 
comprise a toolbox to convert (tensorial) systems of partial differential evolution equations to par- 
allelized C or Fortran code. Kranc can be used as a "rapid prototyping" system for physicists or 
mathematicians handling very complicated systems of partial differential equations, but through 
integration into the Cactus computational toolkit we can also produce efficient parallelized produc- 
tion codes. Our work is motivated by the field of numerical relativity, where Kranc is used as a 
research tool by the authors. In this paper we describe the design and implementation of both the 
' Mathematica packages and the resulting code, we discuss some example applications, and provide 

' results on the performance of an example numerical code for the Einstein equations. 
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I. INTRODUCTION 



The numerical solution of partial differential equations (PDEs) is a common and often quite challenging problem 
in physics, applied mathematics and other fields. The Kranc software presented here is designed as a tool for PDE 
^ ' systems where complexity is a principal problem, and manual coding of equations and general equation-dependent 
CO ' code is too error prone or simply not feasible. Our software is released under the GPL (GNU General Public License), 
. and is available at the URL http://numrel.aei.mpg.de/Research/Kranc Currently, Kranc only supports the 
numerical solution of initial value problems. Extension to initial boundary value problems and elliptic boundary value 
problems is planned for the future. 

Our work is motivated by the field of numerical relativity, i.e. the numerical solution of the Einstein equations j^]. 
(•~^ ', An important open research problem is to find an evolution system for the Einstein equations which is well-suited for 
' numerical simulations. The presented software drastically reduces the complexity of the task of comparing different 
formulations of the Einstein equations; these are complicated nonlinear systems which can have dozens of evolution 
I variables and thousands of source terms. 

^jj^' We aim to address the problems of handling complexity and of dealing with tensorial (or other multi-component) 
• • I quantities in a transparent way. To this end, we present a suite of Mathematica packages, termed "Kranc" , to 
. ^ ' convert (optionally tensorial) systems of partial differential evolution equations to parallelized C or Fortran code. Our 
approach is to use Mathematica to perform high level tensor operations and convert the equations into component 
^ form. This form is then converted to C or Fortran 90. See Sec. Ill Al for the rationale behind this choice. For 
^ evolution problems in three spatial dimensions performance is critical. The design decisions that we have made 
(see Sec. IIIIIII|I allow us to automatically generate codes with efficiencies comparable to hand-written code. Kranc 
provides a simple programming model for code generation. Modifying Kranc to generate more efficient code by-for 
example-re-arranging loops, is straightforward. 

In numerical relativity, the use of computer algebra (CA) methods based on symbolic tensor calculus can not only 
save scientists valuable time, but can also help to substantially reduce errors. It can be used not only for manipulating 
the evolution equations, but also for deriving constraint propagation systems and perturbation formalisms. Having 
a code-generation system directly compatible with the CA system is obviously desirable. Generally, we feel that our 
use of CA makes it easier to focus on algorithms, detached from a particular system of equations. Stressing a more 
abstract point of view is not only mathematically more appealing but also increases flexibility, which benefits scientific 
productivity. 

Our paper is organized as follows: In Sec. ^ we describe the computational infrastructure that we use. Since our 
software only targets the issue of complexity of the equation systems and of processing tensorial quantities, we rely on 
external software to provide a general computer algebra environment and a framework for parallelization. The design 
of the generated codes is described in Sec. IIIII while the design of the Kranc software itself is described in Sec. IIVI 
In Sec. El we describe the use of Kranc to generate codes for the Klein-Gordon equation, the Maxwell equations, 
and a simple evolution system for general relativity, the ADM equations. We consider our software both as a "rapid 
prototyping" system and as a system to create efficient production code, we therefore also include some results on 
performance. 
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II. CHOICE OF INFRASTRUCTURE, AND RELATED DESIGN DECISIONS 

A. The role of computer algebra 

When designing a numerical code to solve tensorial equations, two principal alternatives come to mind. A first 
possibility is to implement high-level tensor operations in a language like CH — h or Fortran 90 using tensorial operations 
like contractions which operate on tensor objects. An interesting example of this approach has been implemented by 
Schnetter . The second option is to use computer algebra to expand all tensorial expressions into components, and 
then use only standard arithmetics in the numerical code. Both approaches have their advantages. While the first 
approach avoids mixing of technologies, the second approach uses specialized technologies to deal with parts of the 
problem. 

We decided to adopt the second strategy for the following reasons. Most importantly, we needed a framework that 
extends beyond code generation, and can also be used as a tool for the analytical part of our work. This requirement 
naturally leads to the use of a computer algebra system. Including the transition from tensors to components and 
the process of code generation in the computer algebra system is relatively easy to implement, and greatly simplifies 
the structure of the resulting numerical code, which is written purely at the level of components. At least part of this 
task could easily have been implemented in a language like PERL, which would have decreased the dependence on 
a commercial computer algebra system (we chose Mathematica, see Sec. Ill B)l . However, this would have demanded 
knowledge of (at least) three different languages from many users and developers: Mathematica, PERL and C/Fortran, 
which would have significantly steepened the learning curve. 

Our design, where the numerical code is not based on tensorial data structures or operations, simplifies performance 
optimization as the code is readable and transparent down to the level of standard arithmetics. The resulting code 
can in principle be implemented in any suitable language, and we provide the option of generating either Fortran 90 
or C (all examples below will be C). Furthermore, in our approach the full power of the chosen computer algebra 
system can be utilized in the process of code generation, e.g. to perform optimizations which rely on the structure of 
the equations. Generally, manipulations on the computer algebra level make it much easier to perform experiments 
involving significant rewrites of the numerical code. 

B. Computer algebra software 

All of the tasks performed by our software which are related to the tensorial nature of equation systems can in 
principle be accomplished by component calculations. These can be carried out quite conveniently and efficiently by 
computer algebra systems like GRTensorll which allows expressions to be entered in an abstract notation and 
yields results in component form. This is often what one wants, but when considering deriving and analysing different 
systems of equations, as is often necessary in numerical relativity research, it is clear that more is needed. Apart 
from the fact that this approach may result in very large calculations requiring significant time and memory, this 
method is unwieldy and not very intuitive. Rather one would like to keep an abstract notation as long as possible, 
and in particular get results in this notation. We have decided to use abstract index tensor calculus; i.e. tensors are 
manipulated symbolically according to rules based on an index notation (see e.g. Wald Q, Sec. 2.4, and Penrose and 
Rindler 0). 

We chose Mathematica as our base CA system for several reasons. The authors were already familiar with its use, 
the language has proven remarkably stable, and Mathematica provides excellent support for pattern matching in the 
core language. Also, with the freely available Ricci j^, xTensor 0] and the commercial MathTensor there are 
already three add-on packages for abstract index tensor manipulation available. xTensor has only been released very 
recently and was not available during the main development phase of the present project. Note that pattern matching 
seems quite essential for tensor manipulations, e.g. Tab and Ted are not the same expression but still are equivalent 
mathematical objects, which can easily be identified with pattern matching techniques. 

Code generation with Kranc can be performed using either MathTensor or the Tensor Tools package, which was 
developed by I. Hinder during the course of this work as a non-commercial alternative. MathTensor has also been 
made the basis of our extensive analytical calculations, but since this article is restricted to the code generation part 
of our project, we will discuss only the usage with TensorTools, which is included in the distribution of Kranc. Both 
MathTensor and TensorTools represent tensors as plain functions of indices: h[la, lb] hat, CD [Metricg [la, 
lb] , Ic] — > 0. This straightforward syntax is less error-prone than Riccis corresponding h [L [a] , L [b] ] or li [L [a] , 
L [b] ] [L [c] ] for a covariant derivative. 

TensorTools uses approximately the same syntax as MathTensor, and will accept reasonably straightforward Math- 
Tensor expressions involving partial derivatives and tensors. This means that tensorial analytical work can be done 
using MathTensor where TensorTools is inadequate, and the results can be fed into TensorTools and hence Kranc. 
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C. Parallel computing infrastructure 

Aiming at "production quality" evolution codes in three spatial dimensions, we regarded parallelization as a vital 
feature. The primary requirements for the parallelization strategy were a short development time and portability - our 
code should run both on inexpensive small clusters, e.g. commodity Linux workstations, but should also scale up to very 
large runs performed at supercomputer centres. We therefore did not want to rely on shared memory or proprietary 
solutions. The current de facto standard for distributed (scalable) computing is MPI (message passing interface) 
0. Several software packages are available which introduce a software layer between the application programmer 
and MPI, and thus significantly reduce the effort to write parallel applications. Two examples we considered as the 
basis for our work are Cactus fl^ and PETSc IT]. PETSc is a general purpose tool developed at the Argonne 
National Laboratory as parallel framework for numerical computations, and Cactus is an open source problem solving 
environment originating in the numerical relativity community, originally developed at the Max Planck Institute for 
Gravitational Physics. While PETSc offers more support for numerical algorithms, in particular for para llel elliptic 
solvers, Cactus already contains some general numerical relativity functionality like horizon finders |l2l IT^ and is 
geared toward flexibility and collaborative projects (it also allows the use of PETSc as a library). We finally decided to 
base our code on Cactus, also because two of the authors currently work at the Max Planck Institute for Gravitational 
Physics, and thus have a chance of direct collaboration with the Cactus development team. 

Cactus mainly targets the issues of parallelization, modularization and portability. The name Cactus derives from 
the design of a central core (or "flesh") which connects to application modules (or "thorns") through an extensible 
interface. Thorns can implement problem specific code on any level from concrete physics applications to low-level 
infrastructure. A set of thorns comprising the "Cactus computational toolkit" provide infrastructure such as parallel 
I/O, data distribution and checkpointing. Current data distribution implementations are based on the principles of 
domain decomposition and message passing using the MPI standard 0. This choice yields an open and reasonably 
documented infrastructure, parallelization, a variety of I/O methods and allows easy interfacing with a growing 
community writing numerical relativity Cactus applications. Modifications to interface with other systems with 
capabilities similar to those of Cactus or with standalone codes should be straightforward. 

We now describe those aspects of Cactus which are important for the remainder of this paper. Cactus does not 
have the structure of a library which provides a set of functions that can be called by the user application. Instead, all 
Cactus thorns are compiled into libraries, and management tasks such as code execution and allocation of distributed 
memory are handled by Cactus and steered via configuration files. The end user supplies a "parameter file" which 
specifies all user-controllable aspects of the run. 

The basic module structure within Cactus is called a "thorn" . All user-supplied code is organized into thorns, which 
communicate with each other via calls to the flesh API (application programmer interface) or APIs of other thorns. 
The integration of a thorn into the flesh or with other thorns is specified in configuration files which are parsed at 
compile time. Thorns are logically grouped into arrangements for organisational purposes. The arrangements live in 
the arrangements subdirectory of the main Cactus directory. Inside an arrangement directory there are directories 
for each thorn belonging to the arrangement. One current shortcoming of Cactus is that this arrangement directory 
structure can only be one level deep. Arrangement and thorn names must be (case independently) unique, must 
start with a letter, and can only contain letters, numbers or underscores. In addition, thorn names must contain 27 
characters or less. 

A key concept for a thorn is the implementation, which defines a group of variables and parameters which are used 
to implement some well-defined functionality. Relationships among thorns are all based upon relationships among the 
implementations they provide. Different thorns providing the same implementation can be compiled into the same 
binary; the decision for which thorn is executed ("active") can be made at run time. 

A thorn consists of a subdirectory of an arrangement containing at least three administrative files, illustrated 
here by examples taken from a Kranc-generated arrangement to solve the Maxwell equations (see the example in 
Sec. lVB|) . The complete code is provided with the Kranc distribution. These files are text files written in the Cactus 
Configuration Language (CCL). The CCL syntax is case independent, and the character indicates that the rest 
of the line is a comment. 

• interface.ccl 

This defines the implementation the thorn provides, and the variables the thorn needs, along with their visibility 
to other implementations. 

# excerpt from Example code EM/EMBase/interf ace . ccl 
implements: EMBase 
inherits: Grid 
public : 
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CCTK_REAL El type=GF timelevels=2 # the electric field 
{ 

Ell, E12, E13 # 3 components of the vector El 

> "El" 

CCTK_REAL Elrhs type=GF timelevels=l # RHS for El-evolution eq. 
{ 

Ellrhs, E12rhs, ElSrhs 
} "Elrhs" 

# analogous for the magnetic field B ... 

• parsun. ccl 

This defines the thorn parameters along with their visibihty to other implementations. The parameters are 
numbers, strings and switches which are read at runtime and control the thorn behaviour. 

# excerpt from Example code EM/EMsetlD/parcim. ccl 
restricted: 

CCTK_REAL sigma "sigma" # a floating point parameter called sigma 

{ 

*:* :: "no restrictions" # all values are allowed 
} # default value is 

• schedule . ccl 

By default no routine of a thorn will be run. Those routines which should be run by the Cactus infrastructure are 
registered in the file schedule, ccl along with directives for appropriate memory allocation and interprocessor 
communication for grid variables. Functions are scheduled in scheduling bins - slots in the main loop such as 
BASEGRID (for setting up coordinates), EVOL (the evolution step) or ANALYSIS (for analysing data). 

# excerpt from Example code EM/EMBase/schedule . ccl 

schedule EMTTBase_RegisterSymmetries at BASEGRID # BASEGRID is a scheduling bin 
{ 

LANG: C 
} "register symmetries" 

III. DESIGN OF KRANC EVOLUTION CODES 
A. Method of Lines 

Aiming both at simplicity and flexibility, we decided to restrict time evolution algorithms to those compatible with 
the method of lines (see e.g. ^3)- Here, the discretization of a system of PDEs describing an evolution problem 
is split into a separate discretization of "space" and "time" . Discretizing only the spatial derivatives (e.g. by finite 
difference, spectral or finite element methods), while keeping the time derivatives continuous yields a coupled set of 
ordinary differential equations (ODEs) in time, with one ODE per spatial grid point. These ODEs are coupled via 
the spatial discretization. For a stable spatial finite differencing scheme, any stable ODE integrator can then be used 
to time-integrate this large system of ODEs. Further key advantages of the method of lines are that: 

• the stability theory is well developed, and is simplified by splitting the analysis into the analysis of spatial and 
time discretizations; 

• numerical methods for the method of lines are well-developed; and 

• by decoupling the spatial finite differencing from the time integration, the resulting numerical scheme and 
computer code are simplified and can easily be modularized. 

In particular, the method of lines makes it straightforward to combine options for spatial discretization and for time 
integration. For the concrete implementation, we choose to base our code on the MoL 15] method of lines thorn within 
Cactus, developed by Hawke. This code provides a parallel ODE integrator, implementing generic Rungc-Kutta and 
iterative Crank-Nicolson methods. 
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B. Spatial Discretization 

Within Kranc, spatial discretization is implemented in a centralized and flexible way. At the moment, only a small 

number of finite differencing methods arc implemented for spatial discretization, but extension to other methods based 
on a regular grid is straightforward. All finite difference formulas are collected in a single header file, currently via 
C-style macros. 

For an evolution code in three spatial dimensions, which evolves equations as complicated as the Einstein equations, 
performance issues are critical. Is is typical that a code spends most of its time evaluating the right-hand-side 
expressions of the evolution equations. It is thus crucial to simplify these expressions as much as possible (this is 
done by the computer algebra part) and optimize the evaluation of the finite difference formulas. It turns out that 
significant performance gains can be obtained when finite difference expressions (and any sub-expressions that are 
very lengthy or re-used many times) are promoted to local scalar variables that are prc-computed at the start of the 
loop over grid points. Apart from cutting down runtime, this method often drastically reduces compile time when 
heavy optimization is used. Pre-computing is the default behaviour. 

The user of the Kranc system will denote derivatives using an abstract notation when passing equations to the 
Kranc functions; e.g. a first derivative of the variable / in the x-direction would be denoted as Dl [f ] . The complete 
list of abstract operators for first and second derivatives is {Dl , D2, D3, Dll, D22, D33, D21, D31, D32}. 

At compile time, the user has two choices to make: 

• Should the derivatives be precomputed at the start of the grid loop? 

• What sort of finite differencing should be performed? 

In order to disable prccomputation, the user should set the preprocessor macro NOPRECOMPUTE, typically through the 
compiler option -DNOPRECOMPUTE. The choice of finite differencing is controlled by defining macros: FD_C2 for second 
order centred, FD_C4 for fourth order centred, or FD_CO for zeroth order (a projection of all derivative terms to zero, 
useful for debugging and performance profiling). These macros may be defined by compiling with, e.g. -DFD_C2. This 
corresponds to the default setting. 

The GenericFD thorn defines several preprocessor macros for performing finite differencing. The concrete finite 
differencing macro name is formed by appending the method to the abstract name; for example, Dl_c2[f] denotes 
the first derivative of / in the x-direction, using centred second order differences. At the moment only second and 
fourth order centred methods are implemented. 

In the C or Fortran code which is produced by Kranc, the abstract notation is replaced with a grid function style 
notation; e.g. Dl[f] is replaced by Dl[f ,i,j,k], where {i,j,k} are the loop indices representing the grid points. 
These "derivative operators" are in fact preprocessor macros, which have different definitions depending on whether 
or not precomputation is enabled. When precomputing is switched off, macros such as Dl [f , i, j ,k] are defined as 
Dlgf [f , i , j , k] . When precomputing is switched on, macros such as Dl [f , i , j , k] are defined as Dlf . Here Dlf is 
used as a local variable which is pre-computed at the start of the loop over grid points; e.g. 

Dlf = Dlgf(f,i,j,k); 

The operators with suffix gf are defined in terms of whichever concrete differencing macros have been selected. 

C. Boundary Conditions 

Currently, Kranc only fully supports periodic boundary conditions. This is achieved via the use of "ghost points" 
on the sides of the grid which Cactus automatically populates with data from the opposite side so that finite difference 
operators can be used right up to the last "real" grid point. 

D. Code organisation 

On the user level, the Kranc code generation toolkit consists of Mathcmatica functions that generate an entire code 
module {thorn in the Cactus terminology). Within the Kranc system, thorns are classified into a small number of 
thorn types, and complex codes can be composed from a larger number of thorns which fall into one of these classes. 
The basic tasks of a thorn are: 

• define Cactus grid functions; 
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• assign values to grid functions; 

• set grid function attributes (e.g. symmetries or boundary conditions); and 

• define Cactus parameters; 

Wc define five types of thorn with specific characteristics regarding their functionahty: the base, setter, translator, 
evaluator and MoL thorns. While the base thorn only defines parameters, grid functions and their properties, all the 
other thorns also perform one or more calculations, i.e. they assign new values to one or more grid functions. Here 
the translator, evaluator and MoL thorns are in a sense "added value" versions of the basic setter thorn. 

A calculation is one of the essential concepts of Kranc. The end result of a calculation is to assign new values to 
a grid function in a loop over grid points. In order to simplify lengthy calculations, all calculations are carried out 
in terms of local scalar variables, and the assignment to grid functions is made at the end of the loop. This avoids 
repeated multiplications due to array accesses. Intermediate expressions, which are only needed as local scalars (i.e. of 
which no derivatives are required), are called shorthands in the Kranc context. 

As an example of a calculation, we list a slightly abbreviated (to save space) version of the file MKGMoL_CalcRHS . c 
(in directory Examples/KleinGordon/MKG/MKGMoL/src), which assigns the right hand sides for the MoL evolution of 
the Klein-Gordon equations. 

// several header files are included and some macros get defined 

/* Define macros used in calculations */ // C++ style comments do not appear in original code 
#define INITVALUE (42) 

void MKGMoL_CalcRHS (CCTK_ ARGUMENTS) 
{ 

DECLARE_CCTK_ARGUMENTS 
DECLARE_CCTK_PARAMETERS 

/* Declare the variables used for looping over grid points */ 

int i = INITVALUE; // saime for j, k, index, istart, jstart, kstart, lend, jend, kend 
/* Declare finite differencing variables */ 

CCTK_REAL dx = INITVALUE; // same for dy, dz, dxi, dyi, dzi, hdxi, hdyi, hdzi 
/* Declare shorthands */ // none in this excimple 
/* Declare local copies of grid functions */ 

CCTK_REAL phiL = INITVALUE; // same for phitL, phirhsL, phitrhsL, Dllphi, D22phi, D33phi 

/* Initialize finite differencing variables */ 
dx = CCTK_DELTA_SPACE(0) ; 
dxi = 1 / dx; 

hdxi = 0.5 * dxi; // einalogously for x and y directions 

/* Set up variables used in the grid loop with stencils suitable for finite differencing */ 
istart = cctk_nghostzones [0] ; 

lend = cctk_lsh[0] - cctk_nghostzones [0] ; // analogously for jstart, kstart, jend, kend 

/* Loop over the grid points */ 
for (k = kstart; k < kend; k++) 
{ 

for (j = jstart; j < jend; j++) 
{ 

for (i = istart; i < lend; i++) 
{ 

index = CCTK_GFINDEX3D(cctkGH,i, j ,k) ; // a Cactus macro to compute the 3D array index 

/* Assign local copies of grid functions */ 
phiL = phi [index] ; 
phitL = phit [index] ; 
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/* Precompute derivatives */ 

Dllphi = Dllgf (phi,i, j ,k) ; // analogous for D22phi and D33phi 

/* Calculate grid functions */ 
phirhsL = phitL; 

phitrhsL = DlKphi.i, j ,k) + D22(phi , i , j ,k) + D33(phi,i, j ,k) - 
phiL*SQR(mass) ; // use macro SQR for mass * mass 

/* Copy local copies back to grid functions */ 
phirhs [index] = phirhsL; 
phitrhs [index] = phitrhsL; 

} 

} 

} 

> 

A Cactus arrangement will typically contain one base thorn, which defines the basic quantities of the model to be 
studied, and an assortment of several of the other Kranc thorn types which implement the equations. Functionality 
which is required but not presently covered by the code that is generated by Kranc, e.g. more involved boundary 
conditions, can easily be added as hand coded thorns. There are currently several restrictions as to what Kranc is 
able to generate: only a few boundary conditions arc implemented, and only periodic boundaries are fully tested and 
supported. 

Note that all thorns can define their own parameters, which can then be shared with other thorns. Apart from 
defining parameters and displaying banner information, the Kranc thorns implement the following tasks: 

• Base thorn: 

— Register grid functions with Cactus. Discriminate between grid functions that can be evolved and need 
two sets of storage (timelevels) and those that only need one time level allocated. 

— Register grid function symmetries. The transformation behaviour under coordinate reflections — > —a;* 
needs to be defined in order to use Cactus infrastructure to run in bitant /quadrant /octant mode for data 
which are compatible with such symmetries. Via the Cartoon method (16j this infrastructure can also be 
used to handle axisymmetric situations in Cartesian coordinates while keeping manifest axisymmetry. 

• Setter thorn: 

— Set grid functions to values, e.g. to set initial data or to update auxiliary grid functions. 

Functions are scheduled in the POSTINITIAL and EVOL time bins at appropriate times relative to the 
time stepping functions. 

• MoL evolve thorn: 

— Register grid functions with the MoL thorn. During registration, the grid functions need to be designated 
as either "evolved" or "primitive" , with respect to this thorn. Evolved variables are those that this thorn 
will evolve using MoL. Primitive variables are those that are used by this thorn in the calculation of the 
right hand sides of the evolution equations, but are not evolved themselves. 

— Generate code to apply boundary conditions. 

— Set MoL right hand side variables through a function scheduled in the EVOL bin 

• Translator thorn: 

— This is a specialized version of a setter thorn which translates between Kranc variables and some other set 
of variables. The idea is that initial data may be specified in one set of variables by an external thorn, and 
can be translated by a translator thorn into the variables you wish to evolve using Kranc. The evolved 
variables can then be translated back into the original set, so that external thorns can use them. These 
translations occur in the POSTINITIAL and POSTSTEP bins, respectively. 

This interface is essential for compatibility with many of Cactus' existing numerical relativity thorns, which 
use the ADM variables via the ADMBase thorn. 
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Mathematica 
script written 
by user 



KrancThorns 

CreateBaseThom[] 
CreateEvaluatorThom[] 
CreateMoLThom[] 
Create SetterThom[ ] 
CreateTranslatorThom[ ] 

Provides functions for 
constructing specific 
types of Kranc thorn 



TensorTools 

MakeExplicit[] 
DefineTensor[] 

Performs manipulations on tensor 
expressions and converts them to 
components 



Thorn 

CreateThom[] 

Constructs a Cactus 
thorn given a 
description in a high 
level form 



CodeGen 

Conditional[] 
GridLoop[] 

Contains functions 
for generating 
blocks of program 
code in a high level 
way 



MapLookup 

Lookup[] 

LookupDefault[] 

MapContains[] 

Contains functions for parsing 
BCranc data structures 



FIG. 1: Relationships between Kranc packages: Each block represents a package, with the main functions it provides indicated 
with square brackets. An arrow indicates that one package calls functions from another 



• Evaluator thorn: 

— Define new grid functions. 

— Assign values to new grid functions in routines scheduled in the Cactus ANALYSIS bin. 

When defining grid functions, tensor components are grouped together as Cactus groups. The reflection symmetries 
of grid functions are read off from their tensor properties automatically (through the function calcSymmetry [gf ] in 
file Thorn. m). By default, grid fimctions are treated as scalars, but e.g. a grid function called hll is interpreted as 
the 11-component of a 2-tensor, which defines its parity information. 



IV. KRANC DESIGN 



Kranc is composed of several Mathematica packages (textual .m files); each one performs a distinct function. The 
user only needs to be concerned with calling functions from one package: KrancThorns. This package contains 
functions for creating the different types of Kranc thorn. 

The diagram in Fig. IIVI illustrates the relationships between the Kranc packages KrancThorns, TensorTools, 
CodeGen, Thorn and MapLookup, which arc described in the following subsections. Structuring the code in this 
way, separating out the different logically independent functions, promotes code reuse. For example, none of the 
thorn generation packages need to know anything about tensors, and none of the packages other than CodeGen (in 
principle) need to know anything about the programming language that is being written (C or Fortran). We have 
chosen to define the concept of a "setter thorn" , and "evaluator thorn" etc, but the mechanics of producing a thorn 
implemented in Thorn and CodeGen are completely independent of this decision. 



A. KrancThorns: Constructing the different Kranc thorns 



We have already discussed the different types of Cactus thorn used in a Kranc arrangement (MoL, setter, base, 
translator and evaluator). The KrancThorns package provides functions to create thorns of these types given high 
level descriptions. These Create*Thorn functions are the ones directly called by users; see the appendix 1^ for 
detailed argument descriptions. Internally the KrancThorns package uses the Thorn package to create the Cactus 
thorns. 
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1. Types of arguments 

Mathematica allows two types of arguments to be passed to a function, "positional arguments" and "named 
arguments" (referred to in the Mathematica book as "optional arguments" . It is possible for some named arguments 
to be omitted from a function call; in this case a suitable default will be chosen. Positional arguments are useful 
when there are few arguments to a function, and their meaning is clear in the calling context. Named arguments are 
preferred when there are many arguments, as the argument names are given explicitly in the calling context. 

For each type of Kranc thorn, there is a function to create it (Create*Thorn). There is a certain set of named 
arguments ("Common named arguments") which can be passed to any of these functions (e.g. the name of the Thorn to 
create, where to create it, etc). Then, for each type of thorn, there is a specific set of named arguments specifically for 
that thorn type. All of the functions accept some positional arguments as well. See the appendix for descriptions 
of the arguments that can be given to each of the KrancThorns functions. 



2. Common data structures 



Kranc is a relatively complex system. It consists of several packages which need to pass data between themselves in 
a structured way. Mathematica does not have the concept of a C++ "class" or a C "structure" , in which collections 
of named objects are grouped together for ease of manipulation. Instead, we have adopted the definition of a "Kranc 
structure" as a list of rules of the form key -> value. We have chosen to use the "rule" symbol "->" for syntactic 
convenience. For example, one might describe a person using a "Person" structure as follows: 

alice = {Name -> "Alice", Age -> 20, Gender -> Female} 

Once a structure has been built up, it can be parsed with the lookup function in the MapLookup package, 
lookup [structure , key] returns the value in structure corresponding to key. For example, lookup [alice , Age] 
would return the number 20. This usage mirrors what is known as an association list (or alist) in LISP style lan- 
guages. Based on this concept we have defined a number of data structures which will be used to describe the thorns 
to construct. Each of these data structures is introduced below, and full details are given in the appendix. 



GroupDefinition 



A GroupDefinition structure lists the grid functions that are members of a specific Cactus group. A list of such 
structures should be supplied to all the KrancThorns functions so that Kranc can determine which group each grid 
function belongs to. 



Calculation 



Calculation structures are the core of the Kranc system; a detailed specification is given in appendix IjA l|l . The 
idea is that the user provides a list of equations of the form variable — > expression. When the calculation is performed, 
for each point in the grid, expression is evaluated and placed into the grid function variable. Here expression may 
contain partial derivatives of grid functions written using the following notation: 

dg r -, dg 

— Di [g] — ^— Dij [g] 

dx^ dx'^dx^ 

where i and j must be given explicitly as integers 1, 2 or 3 representing x, y or z. 

You may specify intermediate (non-grid) variables called "shorthands" which can be used as variable for precom- 
puting quantities which you will use later in the calculation. To identify these variables as shorthands, they should 
be listed in a "Shorthands" entry of the Calculation. 

The arrangement of the terms in the equations can have a marked effect on both compile time and run time. It is 
often helpful to tell Mathematica to collect the coefhcients of certain types of term, rather than expanding out entire 
expressions. To this end, the user can include a "CoUectList" entry in a calculation; this is a list of variables whose 
coefficients should be collected. 

There is the facility for performing multiple loops in a single calculation structure; this can be used to set a grid 
function in one loop, then evaluate derivatives of it in a later loop. For this reason, the equations are given as a list 
of lists of equations. Explicit synchronization of ghost zones is performed after each loop for those groups whose grid 
functions have been set. 

Note that the system is not designed to allow the same grid function to be set more than once in a single loop of a 
calculation. 
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GroupCalculation 

A GroupCalculation structure associates a group name with a Calculation which is used to update the grid 
functions in that group. This is used when creating evaluator thorns, where the calculations are triggered by requests 
for output for specific groups. 

B. TensorTools 

1. Overview 

The TensorTools package was written specifically for the Kranc system, though it is in no way tied to it. We needed 
to perform certain operations on tensorial quantities, and there was no free software available which met our needs. 
TensorTools has the following features: 

• Expand covariant derivatives in terms of partial derivatives and Christoffel symbols (multiple covariant deriva- 
tives can be defined) 

• Expand Lie derivatives in terms of partial derivatives 

• Automatic relabelling of dummy indices to avoid conflicts 

• Convert abstract tensor expressions into components 

Before using any TensorTools functions, the TensorTools package must be loaded. It is self-contained, with no 
dependencies on any other packages. 

Get ["TensorTools .m"] ; 

2. Representation of tensor quantities 

Tensorial expressions arc entered in mostly the same syntax as is used by MathTensor. An abstract tensor consists 
of a kernel and an arbitrary number of abstract indices, each of which can be upper or lower. Abstract indices are 
alphabetical characters (a-z, A-Z) prefixed with either an 1 or a u depending on whether the index is considered to 
be lower or upper. The tensor is written using square brackets as 

kernel [ indices separated by commas ] 

For example, Tj" would be written as T[la,ub]. There is no automatic index raising or lowering with any metric. 
Before using a tensor, it must be registered with the TensorTools package using the Def ineTensor function: 

Def ineTensor [T] 

Entering a tensorial expression causes it to be rendered in an easy to read form: 
In := T[la,lb] 
Out = Tab 

Internally, tensors are represented as Tensor {.kernel. Tensor Index {.label, type} , . . . ] where label is the alphabetical 
index, and type is either "u" or "1" depending on the position of the index. This representation helps in pattern 
matching, and allows TensorTools to identiiy whether a certain object is a tensor or not. 

3. Expansion of tensorial expressions into components 
The function MakeExplicit converts an expression containing abstract tensors into a list of component expressions: 
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In := MakeExplicit[T[la, lb]g[ub, uc]] 

Out = { gll Til + g21 T12 + g31 T13, gl2 Til + g22 T12 + g32 T13, 
gl3 Til + g23 T12 + g33 T13, gll T21 + g21 T22 + g31 T23, 
gl2 T21 + g22 T22 + g32 T23, gl3 T21 + g23 T22 + gSS T23, 
gll T31 + g21 T32 + g31 T33, gl2 T31 + g22 T32 + g32 T33, 
gl3 T31 + g23 T32 + g33 T33} 

Note here that upper and lower indices are not distinguished between in the component form. Tensor Tools was written 
mainly for automated code generation rather than symbolic manipulation; we suggest using different kernels for the 
different forms if this is a problem. 

4- Covariant derivatives 

TensorTools allows the user to define as many covariant derivatives (and associated Christoffel symbols) as are 
necessary. The following defines a covariant derivative operator CD with Christoffel symbol H: 

Def ineConnection [CD , H] 

The function CDtoPD is used to replace covariant derivatives with partial derivatives in any expression: 

In : = CDt oPD [CD [V [ua] , lb] ] 

Out = y% 6 + H\^V'' 

The function MeikeExplicit will automatically do this before converting expressions into components. 

The fundamental operation is to convert a first covariant derivative of an abstract tensor into a partial derivative 
plus some extra terms involving Christoffel symbols. For convenience, TensorTools understands that covariant deriva- 
tives are linear and satisfy the Leibniz property. Also, higher order covariant derivatives are converted to repeated 
application of a first order derivative. 

We define a number of rules to perform the operations in the reduction to the desired form. In the following, x and 
y represent expressions which may or may not contain tensorial indices. 

• Replace any high order covariant derivatives with repeated application of a first order covariant derivative. This 
ensures that we only need to know how to evaluate a first derivative. 

• Replace the covariant derivative of a product using the Leibniz rule: 

Vo(a:y) {V ax)y + x{V aV) 

• Replace the covariant derivative of a sum using the linearity property: 

Va(a; + y)^ VaX + VaV 

• Replace the covariant derivative of an arbitrary expression containing tensorial indices with its expansion in 
terms of a partial derivative and Christoffel symbols, one for each index in the expression: e.g. 

5. Lie derivatives 

The Lie derivative of an expression x with respect to a vector V is written 

Lie[x,V] 

where V has been registered using Def ineTensor and is written without indices. The function LieToPD is used to 
replace Lie derivatives with partial derivatives: 

In := LieToPD [Lie [T[ua, lb] , V]] 

Out = T\^V^ + T%V^^,-T%V- 

Lie derivatives of products and sums are supported. The function MakeExplicit will automatically perform this 
replacement before converting expressions into components. 
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6. Automatic dummy index manipulation 

When two expressions both containing a dummy index b are multiphed together, one dummy index is relabelled so 
as not to conflict with any other index in the resulting expression: 

In := (T[la, lb]g[ub, uc])v[ub, Id, lb] 

Out = Tabg^'V^a^ 

This requires that every multiplication be checked for tcnsorial operands. This can be a performance problem, so the 
feature can be enabled and disabled with SetEnhancedTimes [True] and SetEntLancedTimes [False] . It is enabled 
by default. 

C. CodeGen: High level code generation 

During the development of the Kranc system, we explored two different approaches to generating Cactus files using 
Mathematica as a programming language. Initially, we used a very straightforward system whereby C statements 
were included almost verbatim in the Mathematica script and output directly to the thorn source file. This approach 
has two main deficiencies: 

• The same block of text might be used in several places in the code. When a bug is fixed in one place, it must 
be fixed in all. 

• It is not easy to alter the language that is produced-for example, it is difficult to output both C and Fortran. 

• The syntax in the Mathematica source file is ugly, with lots of string concatenation, making it difficult to read 
and edit. 

To address the first problem, we decided to create Mathematica functions to represent each block of text. This 
also allows the block to be customized by giving the function arguments. By making this abstraction, it became very 
easy to change between outputting C and Fortran, and the resulting Mathematica scripts are more elegant. We call 
this system "CodeGen" . The effect of migrating our old code QJl to the CodeGen system was that the Mathematica 
source files became dramatically shorter and easier to maintain. 

Fundamental to the CodeGen system is the notion of a CodeGen "block" ; this can be either a string or a list of 
CodeGen blocks (this definition is recursive). All the CodeGen functions return CodeGen blocks, and the lists are all 
flattened and the strings concatenated when the final source file is generated. This is because it is syntactically easier 
in the Mathematica source file to write a sequence of statements as a list than to concatenate strings. 

Many programming constructs are naturally block-structured; for example, C f or loops need braces after the block 
of code to loop over. For this reason, it was decided that CodeGen functions could take as arguments any blocks of 
code which needed to be inserted on the inside of such a structure. 

By default, Kranc generates C code. To generate Fortran code the command SetSourceLanguage'Tortran" has to 
be issued after loading KrancThorns .m. Note that while Kranc-generated C-code should always compile (if different 
thorns have been created consistently), a Kranc generated Fortran code may not compile in certain cases. The reasons 
are that the Fortran standard limits names to 31 characters, and specifies a maximum number of 19 continuation lines 
per statement (typical compilers have much larger continuation line limits). Both limits can be compiler dependent, 
and are not checked by Kranc. 

D. Thorn: Constructing the most general Cactus thorn 

We have decided to define several different types of Kranc thorn (base, MoL, setter, evaluator, translator) and to 
provide convenient interfaces for generating each of them. However, there is a lot that is common to each of these, 
and it would be wasteful for the underlying implementation to have separate code for generating each type of thorn. 
The Kranc system is designed to be as modular as possible, with as much code re-use as possible. To this end, there 
is a package called "Thorn" which takes as input a high level description of a Cactus thorn, and which creates the 
necessary thorn flies. 

The interface. ccl, schedule. ccl and param.ccl files are first created by calling Createlnterf ace, CreateParam and 
CreateSchedule with high-level descriptions of these flies, which are converted into CodeGen blocks representing the 
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CCL files. These blocks are then passed, along with blocks for each source file, to CreateThorn, which creates the 
final thorn. The Thorn package provides every aspect of thorn generation that is not specific to the Kranc thorns we 
have defined, and could be used for creating any type of thorn. 

V. EXAMPLES 

In order to illustrate how the packages TensorTools.m and KrancThorns.ni can be used, we give three examples. 
The first example is the massive Klein Gordon equation which, although very simple, contains the essential steps for 
automatic thorn generation. It is described in some detail to serve as a reference for analogous steps in the other 
examples. As a second (still trivial) example, we consider the vacuum Maxwell equations. Here the evolution variables 
are components of tensors and therefore the concept of groups and a GroupDef inition structure becomes essential. 
Finally, as a non-trivial example, we consider the Einstein equations in the standard ADM form. The automatic 
generation of Cactus thorns for these equations involves (almost) the full functionality of the packages TensorTools.m 
and KrancThorns.m. 

The following descriptions should be understood as accompanying explanations of the scripts 
Examples/KleinGordon/MKGTT.m, Examples/Maxwell/EMTT.m and Examples/ADM/KrancADMTT .m which actu- 
ally perform the automatic generation of code. The reader therefore is strongly advised to read the text along with 
the example scripts. Detailed instructions for using the examples are provided in the file README. 

A. The Massive Klein Gordon Field 

As a simple example we consider the massive Klein Gordon field equation 

□ (/)- 771^0 = 0, (5.1) 

where □ is the wave operator on flat space, □ = —d^ + S^^didj in standard Cartesian coordinates {t, x*). 

1. Continuum Formulation 

Introducing 11 = dt4> reduces Eq. I|5.1|l to a system of PDEs which are first order in time 

dt^ = n, 

dtn = 5'^d,dj(l)-m^(j). (5.2) 

We consider Eqs. (|5.2|l on a 3-torus T (i.e. on a rectangular domain with periodic boundary conditions), with initial 
data 

0(0, X*) = cos(27ra;/d), n(0, x^) = 0, (5.3) 

where d is the range of the coordinate x, —d/2 < x < d/2. 

As an example of a quantity which is evaluated for analysis purposes, we introduce the total energy E = J^pd^x, 
with energy density p, 

p = n'^ + 5'^d^(l>dj(l> + m^(l)'^. (5.4) 
The energy is conserved in time and can be used to monitor the accuracy of the numerical solution. 

2. Specifications for Kranc 

The script MKGTT.m in the directory Examples/KleinGordon generates Cactus thorns to solve this time evolution 
problem. After loading the packages TensorTools.m and KrancThorns.m we specify Eqs. (I5.2|l and (|5.4|l as rules 
written in Tensor Tools syntax; i.e. 

EvolEqs = {dot [phi] -> phit , 

dot[phit] -> PD [phi, la, lb] KD[ua,ub] - mass~2 phi}, 
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where KD denotes Tensor Tools' Kronecker delta and the variable IT has been renamed to phit in order to avoid conflicts 
with the constant tt in a Fortran code. The sums over dummy indices are expanded automatically when applying the 
function MakeExplicit, 

{dot [phi] -> phit, 
dot [phit] -> Dll[phi] + D22[phi] + D33[phi] - mass~2*phi} 

The numerical implementation of the above problem is as follows. The evolution variables phi, phit, and the 
parameter mass are declared in a base thorn. Initially the grid functions phi , phit are set to the values H5.3f) by a 
setter thorn. The right hand sides of the evolution equations (|5.2|) are specified in the MoL thorn. Finally the energy 
density rho 1)5. 4|l is evaluated by an evaluator thorn. 

The evolution variables phi and phit have to be given in form of groups 

MKGvars = {{"phi", {phi}}, {"phit", {phit}}}; 

as arguments to the functions CreateBaseThorn, CreateMoLThorn, CreateSetterThorn and should also be given to 
CreateEvaluatorThorn. There is no need to introduce any primitives and therefore we set 

MGKprimitives = {}; 

which is given as a positional argument to the function CreateBaseThorn. 
The parameter mass 

MKGRealBasePareuneters = {mass}; 

is declared in the base thorn and is used in the MoL thorn and the evaluator thorn. It is listed in the 
RealBaseParsimeters argument to the functions CreateBaseThorn, CreateMoLThorn and CreateEvaluatorThorn. 

Equations (|5.2|) . (|5.3|) and H5.4|l and the corresponding grid functions are specified in calculations. For the evolution 
equations we define 

MoLCalculation = {Equations -> {{dot [phi] -> phit, 

dot[phit] -> Dll[phi] + ...}}}; 

Note that the equations are given as a list of lists (consisting of one element here), each sublist resulting in a separate 
loop in the code. 

The calculations endensCalculation for the energy density and IDCalculation for the initial data are set up in 
an analogous way. In addition we assign the name Name -> "time_symmetricID" to the function that sets the initial 
data. The function "time_symmetricID" in the Setter thorn should be called only at the beginning, which is achieved 
by SetTime -> "initial_only". 

The function CreateEvaluatorThorn takes two positional arguments, EvaluationDef initions and a group struc- 
ture which should contain the groups to be evaluated by the thorn (as well as all the groups who's variables are 
referred to). EvaluationDef initions specifies the calculation endensCalculation which is responsible for updating 
the variables in the group "endens", 

EvaluationGroups = {{"endens", {rho}}}; 

EvaluationDef initions = {{"endens", endensCalculation}}; 

We also specify a name for the arrangement, 
MKGSystemName = "MKG"; 

Summarizing, the following information is given to the functions Create*Thorn: 

baseThornInf o = CreateBaseThorn [Join [MKGvars, MKGprimitives] , 
Map [First , MKGvars] , 
Map [First, MKGprimitives], 
SystemName -> MKGSystemName, 

RealBaseParameters -> MKGRealParameters] ; 

mainMoLThornlnf o = CreateMoLThorn [MoLCalculation, MKGvars, 
RealBaseParameters -> MKGRealParameters, 
PrimitiveGroups -> Map [First, MKGprimitives], 
SystemName -> MKGSystemName, 

ThornName -> MKGSystemName <> "MoL", 
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DeBug -> True] ; 

setlDThornInf o = CreateSetterThorn[IDcalculation,MKGvars, 
SystemName -> MKGSystemName , 
ThornName -> MKGSystemName <> "setID", 
SetTime -> "initial_only" , 
DeBug -> True] ; 



evaluateThornInf o = CreateEvaluatorThorn[EvaluationDef initions, 
Join [EvaluationGroups , MKGvars] , 
RealBaseParameters -> MKGRealParameters , 
SystemName -> MKGSystemName, 

ThornNaime -> MKGSystemName <> "evalEnergydens" , 

DeBug -> True] ; 

As a last step, a thorn list MKG.th is created by the function CreateThornList, 

thorns = Join [baseThornInf o , mainMoLThornInf o , setlDThornInf o , evaluateThornInf o] ; 
CreateThornList [thorns , SystemName -> MKGSystemNsime] ; 

The parameter file for running the code has to activate all necessary thorns. In this case, all the thorns in the file 
MKG . th are necessary. The syntax is 

ActiveThorns = "Boundary CoordBase ..." 

The parameter mass is specified as 
MKGBase: :mass =1.0 

We refer the reader to the Cactus thorn documentation for details of the many standard parameters which can be 
adjusted, for example to set up the grid in a different way, or to change the MoL time integration scheme. Here we 
just mention the parameters necessary for changing the resolution and the Courant factor (At/ Ax). The number of 
grid points in the x-direction is set by driver: :global_nx. This number includes both the physical grid points plus 
the ghost points, which are specified by pugh: :ghost_size. grid: :dxyz gives the grid spacing (the same for x, y and 
z). cactus: : cctkj.nitial_time and cactus: :cctk_itlast set the initial time and the last time step. The Courant 
factor is specified by time: :dtfac. 



B. The Maxwell Equations 



As a second example we consider the vacuum Maxwell equations on a Minkowski background in standard Cartesian 
coordinates (t,x'), 

a^F^"- = 0, = 0, (5.5) 

where F^^, is the electromagnetic field strength. 



1. Continuum Formulation 



An observer at rest with respect to the above coordinates measures an electric and magnetic field 

E, = F,o, B, = e.'^'Fku (5.6) 

where i, fc, I run from 1 to 3 and e^fe is the alternating symbol with £123 = 1. Eqs. (|5.5|) can be decomposed into the 
evolution equations 

dtE, = e^'^djBk 

dtB, = -e.^'^djEk (5.7) 

and the constraints 

Ce = 6'^d^Ej =0, Cb^ S'^d.Bj = 0. (5.8) 
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If the constraints (|5.8|l are satisfied initially the evolution equations (|5.7|l guarantee that they are satisfied for later 
times as well. It is therefore sufficient to solve (|5.7f) subject to initial data satisfying (|5.8|) . 
Again we take the spatial domain to be a 3-torus T and consider the initial data 

i?i(0,.^)=acos(?!^), 

E2{Q,X') = (T - 1 cos — -(TCOS \, ^' ), 

a a 
53(0, x^) = (1 - a) cos(^) + a cos(^^^^^tl) ) 

E3{0,x') = BiiO,x') = B2i0,x') = 0, (5.9) 

where —d/2 < x,y < d/2. For ct = these data evolve as a wave travelling in x direction, E2{t,x') = -cos(2l^^) = 
—B3{t,x^), with all other components vanishing. For a = I the wave evolves in a diagonal direction in the xy plane. 
The energy density p is given by 

p=^S'^iE,Ej+B,Bj). (5.10) 
and gives rise to the total energy E = pd^x, which again is conserved in time. 



2. Specifications for Kranc 

The procedure for generating Cactus thorns is similar to the previous example, the major difference being that the 
variables are components of tensors and equations (|5.7|l have to be split into equations for the individual components. 

The automatic generation of code is performed by Examples/Maxwell/EMTT .m. We define the tensors El and B 
with Tensor Tools' function Def ineTensor, where the electric field has been renamed in order to avoid a name clash 
with the exponential constant in Mathcmatica. The evolution equations (|5.7|l in TcnsorTools syntax read 

EvolEqs = {{dot [El [la] ] -> (Eps [la,lb,lc] KD[ub,ue] KD[uc,uf] PD [B [If ] , le] ) , 

dot [B [la]] -> -(Eps[la,lb,lc] KD[ub,ue] KD[uc,uf] PD [El [If ] , le] ) » ; 

These equations are split into individual components by applying the fmiction MakeExplicit. 

As before, we generate a base thorn, a MoL thorn, a setter thorn for the initial data and an evaluator thorn for 
evaluating the constraints and energy density. 

The evolution variables are the electric and magnetic field 

EMvars = {{"El", {Ell, E12, E13», {"B", {Bl, B2, B3»>; 

As before, we do not introduce any primitives. The parameter sigma entering the initial data (|5.9|) only has 
to be known by the setter thorn. It is therefore given as the named argument RealParameters to the function 
Great eSetterThorn. 

We define two groups of variables to be set by the evaluator thorn 

EvaluationGroups = {{"Constraints", {CEl, CB}}, 
{"endens", {rho»>; 

and correspondingly two calculations that are responsible for evaluating the constraints and the energy density. 



C. The Einstein Equations in ADM form 



As a last example we consider the Einstein equations in ADM form, which provide a nontrivial example that 
illustrates the advantages of automatic code generation. A mere 200 lines of Mathematica input is needed to generate 
a whole evolution code including the evaluation of the constraints. Compared to the previous example, we make use of 
more features of Tensor Tools, for example defining a covariant derivative and dealing with tensors with symmetries. 
The code generation part illustrates the concept of primitive grid functions and the possibility of generating two 
different setter thorns to compute the same quantity in different ways. Further, it explains how grid functions are 
computed that depend on derivatives of other grid functions. 
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1. Continuum Formulation 



We decompose the four dimensional vacuum Einstein equations G^j/ ~ into a set of evolution equations and 
constraints corresponding to the standard ADM formulation given by York . We will not touch upon issues of well- 
posedness, as these are irrelevant for our current purpose. Space-time is foliated by a family of spacelike hypersurfaces 
Tit, labelled by a parameter t. Introducing coordinates {t,x^) adapted to this foliation, the four dimensional metric 
can be written as 

ds^ = _(q,2 _ p^pi)dt^ _ 2(5dtdx' + h.jdx'dx^, (5.11) 

where hij is the induced metric on St, a is the lapse function and the shift vector and where we set (3i = hijP^ . Lapse 
and shift describe the components of the time flow dt orthogonal and parallel to the hypersurface Y,t, dt — an + P'^di. 
The rate of change of the induced metric hij along the normal direction n defines the extrinsic curvature Kij = ^Cnhij. 

Projections of the vacuum Einstein equations in the form ''''i?^,^ = into the hypersurface Ei, together with the 
definition of the extrinsic curvature, yield evolution equations for hij and Kij, 

dthij = Cphij + 2aKij (5.12) 
dtK.j = £pK,j+D,D,a + 2ah'"'KijK„,,-aKK,,-aR,j. (5.13) 

Projection of G^^ = in the orthogonal direction yields the Hamiltonian and momentum constraints, 

Ch = R-h'^'hP^'KipK^g + K^ ^0 
Cau = h'"'DiK„n- D^K = 0. (5.14) 

Ci3 denotes the Lie derivative in the direction of the shift vector /3* and Di is the covariant derivative compatible with 
the 3-metric hij. The three dimensional scalar curvature R, Christoffel symbols ^jk and the trace of the extrinsic 
curvature K have been introduced as abbreviations, 

K = h'^K,^, 
fjk = h'\dkhij +djhki~dihjk)/'2, 

R = h'^R,j. (5.15) 

Again the constraints (|5.14|) are propagated by the evolution equations H5.12|l and it is sufficient to evolve freely 
(i.e. without imposing the constraints at each time step) initial data satisfying H5.12|l . 
We express the Ricci tensor Rij in two ways, one involving second derivatives of the metric 

i?y = h'^"'{-d,djhhn + djdrnhu + d,dlhmj - dld„^h,j)/2 + /lim^^'(7'mi7^9J - ihjl'^raq), (5.16) 

and the other is the standard representation in terms of Christoffel symbols 

R^l = - d,Yp, + l^^Jl%q - (5.17) 

Note that Eqs. H5.16|l and (|5.17(l lead to different stencil sizes when discretized, which might change the behaviour of 
the numerical evolution code. 

For our choice of gauge, initial data and boundary conditions, we follow the "gauge wave testbed" described in [l^ . 
The spatial domain is again taken to be a 3-torus T and initial data are given by 

/ill 0,a:' = 1-Asin(— , 

a 

KMOy) ^ ^4^^= (5.18) 
'l -^sin(2li^) 

with all other components of the extrinsic curvature and all off-diagonal 3-metric components vanishing initially. The 
remaining diagonal metric components are 1 initially. A is a constant amplitude and d again denotes the range of the 
coordinate x. The shift is set to zero and the lapse a is given by the square root of the determinant of the 3-metric 
a = Vdet/i. 

The analytical solution to this initial value problem is a pure gauge wave travelling in the x direction at the speed of 
light, hii{t, x^) = hii{x — t), Kii(t, a;*) = Kii{x — i), where all other evolution variables maintain their initial values. 
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2. Specifications for Kranc 

The conversion of tcnsorial equations to equations for individual components using the package Tensor Tools clearly 
is more involved than in the previous examples. We start by defining the tensors appearing in Eas. (|5.12|l - (|5.15|) 

Map[Def ineTensor,{li, hlnv, K, beta, Ricci, gamma, CM}] 

The inverse metric /i*-' has been renamed to hinv as it is going to be a shorthand computed from the metric hij. 
Ricci and gamma denote the Ricci tensor and the Christoffel symbols respectively, and beta is the shift. 

Note that since the tensors hij and Kij are symmetric, we only evolve a subset of linearly independent components. 
To this effect, before applying MakeExplicit to expressions involving symmetric tensors, the symmetries are regis- 
tered. Note that, in conformance with the conventions of MathTensor, in our example covariant and contravariant 
indices are treated differently with respect to the ordering of coordinate indices. The "12" component of a covariant 
symmetric 2-tensor automatically is converted to the "21" component, whereas for contravariant symmetric ten- 
sors the index ordering is the opposite. The corresponding commands for symmetric (co/contravariant) tensors are 
AssertSymmetricDecreasing and AssertSymmetricIncreasing, 

AssertSymmetricIncreasing [hInv [ua,ub] ] ; 

Map [AssertSymmetricDecreasing, {h[la,lb], K[la,lb], Ricci [la, lb] }] ; 
AssertSymmetricDecreasing[gamma[ua, lb, Ic] , lb, Ic] ; 

The symmetry handling in Tensor Tools is rudimentary, and currently only supports symmetries for pairs of indices in 
any particular tensor. When applying MakeExplicit to one of the above tensors all components are kept thus resulting 
in a list with duplicate entries. These have to be removed by a user defined function which is called MctkeExplicitSet 
here. 

To complete the description of the geometric structure for TensorTools the covariant derivative together with the 
corresponding Christoffel symbol is defined, 

Def ineConnection[CD, gamma] ; 

Computing the determinant of the spatial metric h by the Mathematica component expression hDet = 
Det [MatrixQf Components [h [la, lb] ] ] the inverse metric is given by a series of shorthands 

deth -> hDet, 

invdeth -> 1 / deth, 

hlnv[ua,ub] -> invdeth hDet Matrixlnverse [h [ua,ub] ] 

Here, deth is computed as a shorthand at run time to avoid it being calculated for each component of hlnv [ua,ub] . 
hDet Matrixlnverse [h [ua,ub] ] is the component inverse of h with the common determinant factor cancelled. 

The quantities appearing in Eqs. H5.12|l and (|5.14|l split in the following way into evolution variables, primitives 
and shorthands: Evolution variables are the 3-metric and the extrinsic curvature 

ADMvars = {{"h", {hll, h21, h31, h22, h32, h33», {"K",{K11, K21, K31,K22, K32, K33}}}}; 

All other quantities appearing in (|5.12|) and (|5.14|) have to be either primitive grid functions, which are set in a setter 
thorn, or shorthands defined inside the loop in which they are used. The decision about which quantities are primitives 
is based on the following considerations. Quantities that appear in differentiated form in the evolution equations, 
e.g. lapse and shift, have to be grid functions. Quantities which can be computed in various different ways, e.g. the 
Ricci tensor, can be set in different setter thorns which then can be compiled into the same executable, where one 
of them is "activated" at run time. Quantities that are needed for both the evolution equations and the constraints, 
e.g. the Christoffel symbols, can be set in a setter thorn in order to reduce computational costs. Furthermore all 
quantities whose output is required, e.g. the trace of the extrinsic curvature as specified in the tests in have to 
be grid functions. All of these quantities should be declared as primitives. The complete list of primitives for our 
example therefore is 

ADMprimitives = {{"alpha", {alpha}}, {"beta" , {betal , ...}}, {"Ricci" , {Riccill , ...}}, 
{"gamma" ,{gammalll, ...}}, {"trK" , {trK}}} 

The scalar curvature is computed as a shorthand from the Ricci tensor. 

Eqs. H5.12|) - H5.17|l in MathTensor syntax are read from the file ADMEqs.m and expanded into components with 
MakeExplicit. In order to increase performance once the functionality of TensorTools is not needed anymore we 
switch off the special treatment of products by setting 
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SetEnhcincedTimes [False] 

As a last preparational step we define the components of the inverse metric and the lapse in a CollectList which 
is going to be used to simplify the equations at the level of components, essentially by applying the Mathematica 
function Collect [rhs, CollectList, Simplify] to the right hand sides of assignments. 

The generation of the base thorn, the MoL thorn, the various setter thorns and the evaluator thorn run along the 
same lines as in the previous examples with the following additional features. All primitive grid functions have to be 
set at every time step, so the argument SetTime for all setter thorns therefore has to be "initial_and_poststep". 

The Ricci tensor in the form Eq. (|5.17|l requires taking derivatives of the Christoffel symbols, which therefore have 
to be known on the whole grid before the Ricci tensor is computed. This is done in the code by using two loops, 
where the first loop is used to compute the Christoffel symbols. At the end of the loop the grid functions that store 
the Christoffel symbols are synchronized such that the ghost zones are filled in with data from the opposite side of 
the grid. The Ricci tensor can then be computed in the second loop. In order to generate two separate loops, the 
calculation passed as an argument to the function CreateSetterThorn has to specify the equations in two separate 
(sub)lists, i.e. 

Equations -> {Join [InvMetricEqs , Christof f elEqs] , RicciStandardEqs}; 

Another major difference to the previous examples is the way the initial data is set. We use the existing thorn Exact 
in the arrangement AElThorns to compute the metric and extrinsic curvature as given in Eq. (|5.18|1 . The only thing 
left to do is to convert the initial data from the variables used by the Exact thorn, {gxx, . . . , kxx, ....}, defined 
by the CactusEinstein arrangement, to our evolution variables and correct for different sign conventions in the def- 
inition of the extrinsic curvature, i.e. hll == gxx, . . . , Kll == - kxx, .... This is done in the Translator thorn, 
which implements functions to translate either way: from the CactusEinstein variables to the evolution variables at 
the POSTINITIAL Cactus time bin, and from the evolution variables to the CactusEinstein variables after every time 
step. The latter translation is performed to allow use of existing analysis tools, e.g. in CactusEinstein. The trans- 
lator functions are defined by the Kranc calculations Translator InCalculation and TranslatorOutCalculation 
respectively. 

The thorn list is generated as in the other examples. This time, however, it has to be sup- 
plemented by the thorns CactusEinstein/ADMCoupling, AElThorns/Exact , CactusEinstein/CoordGauge and 
CactusEinstein/StaticConf ormal which are needed by the Exact thorn to set up the initial data Eqs. (|5.18|) . 

All the thorns contained in the thorn list - including both Setter thorns for the Ricci tensor - are compiled into 
one executable. The parameter file specifies the way the Ricci tensor is computed by activating the relevant thorn. 
The specification of the parameters determining the initial data can be found in the documentation of the thorn 
AEIThorns/Exact. 



3. Performance 

We include some timing results for simulations to demonstrate that our computer-generated code can achieve a 
performance that is comparable to carefully optimized handwritten code. We do not attempt a systematic study of 
performance issues, and restrict ourselves to a single example based on the "gauge wave" data as defined in Eq. H5.18|l . 
Our test run is a 3D version of the gauge wave spacetime defined with a grid size of 50"^ (plus one ghost zone); we 
evolve for 100 time steps. 

We compare the performance of the ADM code presented above to a code that has been developed as a production 
code by the AEI astrophysical numerical relativity group. This code is publicly available from the Cactus CVS server 
as the CactusEinstein arrangement [23|. The Kranc and CactusEinstein ADM codes use the same numerical methods 
(3 step Iterative Crank Nicolson and second order centred spatial finite differencing). As a physics production code, 
CactusEinstein ADM is significantly more complex, but the evolution algorithm is essentially equivalent. The numbers 
given below should only be taken as a rough guide. In contrast to Kranc ADM, which uses computer generated code 
and a generic infrastructure for spatial differencing and time evolution, CactusEinstein ADM is hand coded and does 
not support general numerical schemes (although different time update schemes have been implemented as well). 

In the tables IV C 31 and IV C 31 compile time and user runtime have been timed with Linux's time utility and show 
the total number of CPU-seconds that the process spent in user mode. The results of the cactus timing column have 
been obtained with the internal timing infrastructure of Cactus. 

Tests have been performed on two machines, a Dell D600 Laptop running Redhat Linux 9.0, with a 1.6 GHz Intel 
Pentium M CPU and 1 GByte of RAM; and a Dell precision 450 workstation running Redhat Linux 7.2, with a dual 
Intel Pentium Xeon CPU running at 1.7 GHz and again 1 GByte of RAM. All runs have been performed as single 
processor runs. 
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code version 


compile time 


user runtime 


cactus timing 


CactusEinstein ADM, Intel 8 


184 


145.2 


145.6 


CactusEinstein ADM, VAST f90 & gcc 


154 


155.0 


155.2 


Kranc ADM, C, Intel 8 


119 


154.9 


155.4 


Kranc ADM, C, VAST £90 & gcc 


89 


166.8 


167.4 


Kranc ADM, F90, Intel 8 


129 


159.7 


159.5 


Kranc ADM, F90, VAST f90 & gcc 


92 


168.3 


168.8 



TABLE I: Timing results from a Dell D600 Laptop running Redhat Linux 9.0, with a 1.6 GHz Intel Pentium M CPU and 
1GByte of RAM, comparing results obtained with the Intel 8 compilers vs. a combination of Pacific Sierra VAST F90 and gcc 
3.2.2 compilers. Optimization options for the Intel compilers are -03 -xN -ip, and -03 for VAST and gcc. 



code version 


compile time 


user runtime 


cactus timing 


CactusEinstein ADM, Intel 7.1 


303 


188.1 


189.5 


Kranc ADM, C, Intel 7.1 


186 


131.3 


131.1 


Kranc ADM, F90, Intel 7.1 


194 


115.3 


115.7 



TABLE II: Timing results from a Dell Precision 450 workstation running Redhat Linux 7.2, with a 1.7 GHz Intel Pentium 
Xeon CPU and 1GByte of RAM, binaries have been built with the Intel 7.1 compilers with optimization options -03 -xW -ip. 



Note that with the exception of the result for CactusEinstein ADM compiled with Intel 7.1, the Kranc generated 
code is around 6 - 10 % slower. Possible explanations are that the CactusEinstein ADM code partially uses Fortran 90 
array language which allows further optimizations (e.g. the vectorization feature of the Intel compiler) , and the scheme 
for synchronizing grid functions within Kranc is not yet optimal (synchronizing is not necessary if no derivatives appear 
in assignments, as the ghost zones can be updated locally). We speculate that the CactusEinstein result for the Intel 
7.1 compiler might be caused by the compiler reaching internal limits and giving up on certain kinds of optimization. 
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APPENDIX A: DATA STRUCTURE SPECIFICATIONS 



Here we describe in detail the data structures which are used when calling the KrancThorns functions. 
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1. Calculation 



Key 


Type 


Description 


Equations 

Shorthands (optional) 

Name (optional) 
Before (optional) 

After (optional) 


list of lists 

list of symbols 
string 

list of strings 
list of strings 


{loopl, loop2} - Each loop is a list of rules of 
the form variable -> expression where variable is 
to be set from expression 

Variables which are to be considered as 'short- 
hands' for the purposes of this calculation 

A name for the calculation 

Function names before which the calculation 

should be scheduled. 

Function names after which the calculation should 
be scheduled. 



2. GroupCalculation 

A GroupCalculation structure is a list of two elements; the first is the name (a string) of a Cactus group and the 
second is the Calculation to update the variables in that group. 

3. GroupDefinition 

A GroupDefinition structure is a list of two elements. The first is the name (string) of a Cactus group and the 
second is the list of variables (symbols) belonging to that group. The group name can be prefixed with the name of 
the Cactus implementation that provides the group followed by two colons (e.g. "ADMBase::metric"). If this is not 
done, then the KrancThorns functions will attempt to guess the implementation name, usually using the name of the 
thorn being created. 

APPENDIX B: KRANCTHORNS FUNCTION REFERENCE 

Here we document the arguments which can be specified for the functions CreateBaseThorn, CreateMoLThorn, 
CreatcScttcrThorn, CreateTranslatorThorn and CreateEvaluatorThorn. 

Note that we use Mathematica syntax for function-specific section headers. Underscores denote function arguments, 
and OptArguments stands for optional arguments, also referred to as named arguments below. These are given in 
the form my Function [... , argumentName -> argument Value] . 

1. Common Named Arguments 



The following named arguments can be used in any of the Create*Thorn functions: 



Argument 


Type 


Description 


Default 


SystemName 


string 


A name for the evolution system implemented by 
this arrangement. This will be used for the name 

of the arrangement directory 


"MyPDESystem" 


SystemParentDirectory 


string 


The directory in which to create the arrangement 
directory 


u 


ThornName 


string 


The name to give this thorn 


SystemName + thorn type 


Implementation 


string 


The name of the Cactus implementation that this 
thorn defines 


ThornName 


SystemDescription 


string 


A short description of the system implemented by 
this arrangement 


SystemName 


DeBug 


Boolean 


Whether or not to print debugging information 
during thorn generation 


False 
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2. Arguments relating to parameters 

The following table describes named arguments that can be specified for any of the thorns except CreateBaseThorn. 

CrcatcBascThorn is special because it can be used to define parameters which are inherited by each thorn in the 
arrangement, so the arguments it can be given are slightly different. 



Argument 


Type 


Description 


Default 


RealBaseParameters 


list of strings 


Real parameters used in this thorn but defined in 
the base thorn 


{} 


IntBaseParameters 


list of strings 


Integer parameters used in this thorn but defined 
in the base thorn 


{} 


RealParameters 


list of strings 


Real parameters defined in this thorn 


{} 


IntParameters 


list of strings 


Integer parameters defined in this thorn 


{} 



3. CreateBaseThorn [groups-, evolvedGroupNames_, primitiveGroupNames_, OptArguments ] 



a. Positional arguments 



Argument 


Type 


Description 


groups 

evolvedGroupNames 
primitiveGroupNames 


list of GroupDefinition structures 
list of strings 
list of strings 


Definitions of any groups referred to in the other 

arguments. Can supply extra definitions for other 

groups which will be safely ignored. 

Names of groups containing grid functions which 

will be evolved by MoL in any of the thorns in the 

arrangement. 

Names of groups containing grid functions which 
will be referred to during calculation of the MoL 
right hand sides in any of the thorns in the ar- 
rangement. 



h. Named argum.ents 



Argument 


Type 


Description 


Default 


RealBaseParameters 
IntBaseParameters 


list of strings 
list of strings 


Real parameters defined in this thorn and inher- 
ited by all the thorns in the arrangement 
Integer parameters defined in this thorn and in- 
herited by all the thorns in the arrangement 


{} 
{} 



4. CreateEvaluatorThorn[groupCalculations_, groups., OptArguments ] 

a. Positional arguments 



Argument 


Type 


Description 


groupCalculations 
groups 


list of GroupCalculation structures 
list of GroupDefinition structures 


The GroupCalculations to evaluate in order to set 
the variables in each group 

Definitions for each of the groups referred to in 
this thorn. Can supply extra definitions for other 
groups which will be safely ignored. 
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5. CreateMoLThorn [calculation., groups., OptArguments ] 



a. Positional Arguments 



Argument 


Type 


Description 


calculation 
groups 


Calculation 

list of GroupDefinition structures 


The calculation for setting the right hand side vari- 
ables for MoL. The equations should be of the 
form dot Igf] -> expression for evolution equa- 
tions, and shorthand -> expression for shorthand 
definitions, which can be freely mixed in to the 
list. 

Definitions for each of the groups referred to in 
this thorn. Can supply extra definitions for other 
groups which will be safely ignored. 



b. Named Arguments 



Argument 


Type 


Description 


Default 


PrimitiveGroups 


list of strings 


These arc the groups containing the grid functions 
which are referred to but not evolved by this evo- 
lution thorn 


{} 



6. CreateSetterThorn[calculation_, OptArguments ] 

a. Positional Arguments 



Argument 


Type 


Description 


calculation 
groups 


Calculation 

list of GroupDefinition structures 


The calculation to be performed 
Definitions for each of the groups referred to in 
this thorn. Can supply extra definitions for other 
groups which will be safely ignored. 



b. Named Arguments 



Argument 


Type 


Description 


Default 


SetTime (optional) 


string 


"initiaLand-poststep" , "initiaLonly" or "post- 
step_only" 


"initiaLand-poststep" 



7. CreateTranslatorThorn[groups_, OptArguments ] 

a. Positional Arguments 



Argument 


Type 


Description 


groups 


list of GroupDefinition structures 


Definitions for each of the groups referred to in 
this thorn. Can supply extra definitions for other 
groups which will be safely ignored. 
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b. Named Arguments 



Argument 


Type 


Description 


TranslatorlnCalculation 
TranslatorOutCalculation 


Calculation 
Calculation 


The calculation to set the evolved variables from 
some other source 

The calculation to convert the evolved variables 
back into some other set of variables 
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